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1 Introduction 

During this period of performance (January 1, 1992 - December 31, 1992) we have accomplished the inves- 
tigation of the following tasks: 

a. Application of the Non-linear Force-free Field (NLFF) Model to study the active region evolution. 

b. Coronal dynamical responses due to emerging flux including the transition region. 

c. Loss of MHD equilibrium due to foot point motions. 

d. Two-dimensional MHD global coronal model: steady-state streamers. 

The detailed description of these studies are included in the following sections: 

2 Applications of the Non-Linear Force- Free Field (NLFF) Model 
to Study the Active Region Evolution 

Since the publication of the first version of our extrapolation scheme (Wu et al. 1990), called ” Progressive 
Extension Method (PEM)”, further modification was made on the regularization-like technique. This new 
regularization-like technique can be summarized as follows: the expression for averaging given by Eqs. (3.5) 
and (3.6) of the paper by Wu et al. 1990 is modified as follows: 

The regularization-like solution for the magnetic field vector, B x 3 becomes 

Si j = (1 — + 7&i>j ( 1 ) 

where 5^ is the same as Eq. (3.5) from the Wu et al. (1990) paper and 

7 {z) = To exp[a(z/z max )] (2) 

with t„ being arbitrary and to be determined by a trial and error method until the solution converges and 

a = ln{y{z ine ix)/to) - (^) 

By setting our goals on t( z maa: ) and 70 , a is determined. 

Using this modified technique, we performed a number of significant tests by using a complicated an- 
alytical nonlinear force- free solution (Low and Lou, 1990) as the input to the numerical model. Figure 1 
shows a top view of the field lines generated by: (a) an analytical solution (Low and Lou, 1990), and (b) a 
numerical solution obtained from extrapolation. Figure 2 shows three-dimensional field lines corresponding 
to the solutions given in Figure 1. These results clearly demonstrate that the PEM (Wu et al. 
1990) is a reasonable method to obtain nonlinear force-free fields by using vector field data in 
the photosphere. The accuracy of the extrapolation can be shown as a function of height. It is recognized 
from these results that the height of the extrapolation is about one tenth of the horizontal boundary which 
is typically about 30,000 - 50,000 km. 

We then applied this NLFF model to extrapolate the magnetic field configuration using the measured 
photospheric vector field at Beijing Observatory on 1989 March 9 and 11 to demonstrate the capability of this 
model. In order to show the validity of this NLFF model, we also show the computed potential field using 
the present algorithm to compare with the Schmidt’s potential field model. This result is shown in Figure 

3 (for 1989 March 10 at 0600 UT) which indicates that the results obtained from the present algorithm are 
identical to the classical Schmidt’s model. 

Figure 4 shows the nonlinear force-free extrapolation for 1989 March 10 (0600 UT) and 11 (0226 UT) 
respectively. From these results, we notice that, (i) the differences between the potential field and the NLFF 
models are significant. That is, the magnetic field structure represented by the NLFF model is believed to 
more closely resemble the realistic situation as shown in Figure 5. We clearly notice from Figure 5 that the 
loop structures revealed by the NLFF model closely resembled the H-a observations which the potential field 
model cannot show. 
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Figure 1. Comparison between (a) the analytical solution and (b) extrapolated solution using analytical 
data (i.e. data generated by solution given in (a)). 



Figure 2: Three-dimensional representation of field lines obtained by analytical solution (a) and numerical 
model (b) corresponding to the solution given in Figure 1 
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Figure 4: Simulated temporal evolution of NLFF Model solution over a period of about 20 hours. The 
dashed and solid contours are the photospheric vector magnetograph measurements. 





3 Coronal Dynamical Responses Due to Emerging Flux Includ- 
ing the Transition Region 

From the observations, it is inferred that the additional magnetic energy stored in the solar atmosphere very 
likely results from flux emergence from the sub-photospheric layer. It is difficult to model this physics because 
the characteristics of this part of the solar atmosphere (i.e. photosphere/chromosphere, transition region and 

corona) have several order-of-magnitude variations in density, pressure and plasma beta (i.e. (3 — » 

hence little progress has been made up to now. Recently, Wu et al. (1992) developed a self-consistent MHD 
model to address this problem. A more realistic energy equation was used in order to construct this model. 
This energy equation, which includes wave heating, radiative cooling, thermal conduction and joule heating, 
is given as: 

^ + (V ■ V) P + 7 p(V • V + (7 - 1)(V •<?- — - Gmech + Lrad) = 0 (4) 

at cr 

where, p = pressure, V == plasma flow velocity vector, Q — thermal conduction, J = electric current, G fnrc h 
— wave heating, L ra a = radiative cooling, 7 = specific heat ratio and a = electrical conductivity. 

The other equations (i.e. continuity, momentum and the Maxwell equations) are identical to those used 
before (Wu et al. 1983b). The numerical method used for this study is the modified ICED- ALE (Implicit- 
Continuous-Eulerian-Difference-Mesh-Arbitrary-Lagrangian-Eulerian) method given by Wu et al. (1991b). 
The reason for choosing this method is because the Eulerian difference scheme usually used has a limitation 
on the grid size which becomes impractical when the spatial gradient becomes too large. A brief description 
of the numerical results obtained from this model is given as follows: Figure 6 shows the steady-state density, 
temperature, pressure and plasma beta distribution from the surface through the temperature minimum up 
to the lower corona (~ 3500 km) together with grid distribution at 50 s after introduction of emerging 
magnetic flux. These results are equivalent to the empirical model of the Harvard- Smithsonian model. Our 
next objective was to demonstrate that this model can be used to study the momentum and energy transport 
from the solar interior to the corona and thereby, to investigate the physical mechanism of coronal heating 
and solar wind acceleration. To accomplish this objective we introduced an emerging flux at the lower 
boundary and computed the evolutionary state of the plasma properties, velocity and magnetic field. The 
initial state is shown in Figure 6. Figures 7 and 8 show the evolutionary results for velocity, temperature, 
density, pressure, and magnetic field lines at time 500 s and 1000, after introduction of emerging flux. From 
these preliminary results, we found the following features (i.e. a paper is currently in preparation and the 
results were presented at the AAS/SPD Meeting June 6-11, 1992): 

i. The emerging flux leads to the formation of a current sheet at the interface of the old and new magnetic 
fields and to its propagation upward toward the corona; 

ii. The induced plasma flow oscillates vertically at the Brunt- Vaisala frequency with a period of 240 s.; 

iii. Also, we showed that there is no oscillation when the gravity (unrealistically) is ignored,; 

iv. The maximum downward flow (~ 20fcms~ 1 ) occurs in the neighborhood of the legs of the magnetic 
loop which is a typically observed feature. 

4 Loss of MHD Equilibrium Due to Footpoint Motions: A 

Three-Dimensional, Time- Dependent MHD Simulation Model 

Recently Sudan (1991) has demonstrated the phenomena of ’Toss of equilibrium” by using a set of reduced 
incompressible MHD equations. We have used a newly developed three-dimensional, time-dependent, com- 
pressible MHD simulation model (Sun and Wu, 1992) to demonstrate similar results. Parker (1972, 1981), 
in a long series of papers spanning almost two decades, has claimed that the coronal magnetic field, evolving 
in response to smooth continuous photospheric footpoint motions, will not be able to achieve a smooth, 
force-force equilibrium; instead, the field develops tangential discontinuties. It has been thought that these 




-L51UT 



Potential Fields at 15 March 1989 0151 UT 



gure 5. Comparison between the magnetic loop structures seen by H-a hit ere rams and magnetic 
loops derived by the NLFF model as well as the potential field model: (a) the magnetic fields structure 
enve y t le inode at 1089 March 15 at 0315 l' T overlay on the H-a picture, (b) the magnetic fields 
s ructure derived by the potential field model at 1989 March 15 at 0315 UT. Note the regions indicated 
by A and B. It is easy to recognize that the potential field model failed to match the H-a picture 
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Figure 7. Evolution of plasma flow (velocity vectors), (a) at t = 100 sec. (b) at t = 200 sec. (c) at t 
= 300 sec. (d) at t = 400 sec. (e) at t = 500 sec. (f) at t = 600 sec. (g) at t = 700 sec. (h) at t - 800 
sec. (i) at t = 900 sec. (j) at t = 1000 sec. 
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discontinuties could lead to enhanced dissipation in current sheets and, moreover, could provide the energy 
release responsible for coronal heating (i.e. x-ray emissions as observed). 

Our simulation for the loss of MHD equilibrium due to footpoint motion is based on a three-dimensional, 
time-dependent, compressible MHD model using a new NICE (Nimble ImpUcit-Continuous-Eulerian) scheme 
(Sun, 1991; Sun and Wu, 1992). The results of the simulation shown in Figure 9 are explained as follows: 
The sequence of panels, with projected field lines onto the x-z plane, shows that, up to about f-max (this 
” maximum time” is deduced from an analytical solution when the force-free solution reaches the critical 
state), the numerical results closely follow the analytical equilibrium sequence. After (t inax ) the field lines 
start to oscillate with a maximum velocity reaching an amplitude of about 50 kms 1 . This result suggests 
that ”loss of equilibrium” has occurred. Since we are using the ideal MHD equations, no reconnection can 
take place, and an eruption is not to be expected with the infinite amount of overlying flux. This resembles 
the real situation for a small arcade of loops embedded in a large active region. With an estimated coronal 
magnetic Reynolds number of about 1 0 1 3 , and with the numerical results giving no indication of strong 
current concentrations, indeed no significant reconnection can occur on a time-scale of up to 100 Alfven 
scale times (the duration of the simulation). Therefore, the conversion of magnetic energy into MHD wave 
energy may be the only efficient method that is available to shed excess free energy and, thus, for heating 
the corona. Details of this study are described by Marten, Sun and Wu (1992). 


5 Two-Dimensional MHD Global Coronal Model: Steady State 
Streamers 

In this subject, we have made some progress since we first constructed a streamer model (Steinolfson, 
Suess and Wu, 1982). Our motivation to revisit this problem was to extend the outer boundary farther 
away from the Sun (i.e. ~~35 solar radii) and to gain the experience necessary for development of a three- 
dimensional model. Another motivation to develop such a model is the simulation of streamers in support of 
the Ultraviolet Coronagraph and Spectroheliograph (UVCS) and the Large Angle Spectrometric Coronagraph 
(LASCO) instruments on the Solar Heliospheric Observatory (SoHO). These instruments will be able to 
measure the temperature, density, and flow velocity vector in the corona. With model calculations, it will be 
possible, for example, to infer the magnetic vector. The results of the present study are presented by Wang 
et al. ( 1992a, b) and Noci et ai (1992). We shall briefly summarize some of the highlights in the following 
paragraph. 

Figure 10 shows the steady state magnetic field lines for four cases: (a) Dipolar 0 O = 0.5, (b) Quadrupole, 
0 = 0.5, (c) Hexapole, 0 - 0.5, and (d) Dipole, (3 = 0.2. The relaxation times allowed to reach these equilibria 
are: (a) 22.22 hrs., (b) 16.7 hrs., (c) 18.06 hrs., and (d) 19.44 hrs. respectively, where 0 is evaluated at the 
equator of the solar surface. In each plot, four dashed lines are labelled ”A, B, C, or D . These show the 
radial directions used for plotting the variables versus radius in each case. Thus, the 0 — 0.5 quadrupole 
plots will have variables versus radius at the pole (B), at the edge of the polar region (A), through the 
mid-latitude streamer (D), and in the equatorial open region (C), The dashed lines are along the direction 
of the grid. Since there is no grid point either exactly on the equator or exactly at the pole, these lines are 
slightly offset from those positions. The plasma parameters, radial velocity and total magnetic field strength, 
are presented in a pre-print (Wang, Wu, Suess and Poletto, 1992b) which is included in the Appendix. 



Magnetic Field Line Projection on the X-Z Plane for Beta = 0*i0 





Figure 9: Numerical simulation of the evolution of the magnetic field in projection on the x — z plane, which 
is perpendicular to the x — y photospheric plane. 
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Abstract. We describe a two-dimensional time-dependent, numerical, maeneiohydrodynarruc model 
for the determination of the physical properties of coronal streamers fro the top of the transition 
zone ( ft. 5 = 1) to 15 FI?,. Four examples are given: for dipole, quadrupole. and hexapole initial 
field topologies. The computed parameters are density, temperature, velocity, and magnetic held. In 
addition to the properties of the solutions, their accuracy is discussed. We use the model as the basis 
tor a general discussion of the way noundary conditions are soecuied in this and similar simulations. 


I. Introduction 

We present results from a recently-developed numerical model of coronal struc- 
ture. The immediate reasons for a new model were to extend the outer boundary 
farther from the Sun and to gain the experience necessary for development of a 
three-dimensional model. A result of this process has been a close examination of 
the physical details of the solution and how they depend on the way the boundary 
conditions are specified. An immediate application will be the simulation of stream- 
ers in support of the Ultraviolet Coronagraph and Spectroheliograph (UVCS) -and 
the Large Angie Spectrometric Coronagraph (LASCO) on the Solar Heliospheric 
Observatory (SoHO). These instruments will be able to measure the temperature, 
density, and flow vector in the corona. With model calculations, it will be possible, 
for example, to estimate the magnetic field vector. 

Numerical models of coronal structure have been published sporadically, at 
long intervals, over the past twenty years. The first (Pneuman and Kopp, 1971) 
demonstrated the feasibility of such models, treating isothermal flow and arriving 
at the solution by iterating on the electrical currents. However, a more efficient- 
and flexible method is to consider an initial-boundary value problem in which 
the steady state is found holding the boundary conditions constant and allowing 
the solution to relax in time from an essentially arbitrary initial state. Steinolfson, 
Suess, and Wu (1982) applied this later technique to the analysis of a poly tropic 
dipole configuration for a range of plasma 3 (ratio of internal pressure to mag- 
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nenc pressure). Steinoifson ( i 989. 19911 ana Guo et d. ;1992) have used this 
steady-state solution as the basis for studying coronal mass ejections and streamer 
evolution with shear, which can be simulated using a nearly identical numerical 
model. Details of the numerical schemes and results can be found in the referenced 
publications. 

We revisited this problem for the reasons mentioned above. However, we also 
consider that such complex numerical models are rareiy without problems or uncer- 
tainties. When the models are used for analysis of data and for predictions, the only 
reliable validation is to develop an independent model and compare the resuits. 
Even when both (or all) models are fundamentally correct, this process generally 
leads to new or deeper understanding of the problem. In the present case, this is 
precisely what has happened. We have gained a better insight into the physical 
basis of the criteria which should be adopted in specifying boundary conditions. 
The results from this constitute an important part of the present study. 

The physical and numerical simulation is described in Section 2. Section 3 
details numerical models of dipole, quadrupole, and hexapole magnetic fields. 
Section 4 is a discussion of numerical precision of the solution and the boundary' 
conditions, putting the discussion into context with earlier models so far as is 
possible. Section 5 contains our summary and conclusions. 


2. The Physical and Numerical Simulation 


We assume axisymmetric. single fluid, polytropic. time-dependent ideal magneto- 
hydrodynamic flow and perform the calculation in a meridional plane defined by 
the rotational symmetry axis of the magnetic field. The coordinates are (r. 9, o) 
with 6 being the ignorable coordinate. For the magnetic field boundary condition, 
we take the radial field component at the lower boundary to be that given by a 
vacuum dipole, quadrupole. or hexapole potential magnetic field. The flow there- 
fore has reflective symmetry across the equator and the calculation need be done 
in only one quadrant. The equations of motion that describe this flow are: 
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dv T 

dr or r at) r 
The dependent variables are the density. p, the pressure, p, the raaial ana 
meridional velocities. ty and v$, and the radial and meridional magnetic fields, 
B r and Bg. The constants M. z , G, 7 , and p. are the solar mass, gravitational 
constant the poiytropic index and the magnetic permeability. 

These equations are solved in a computational domain extending from the Sun 
(1 Rr v ) to 15 Rq, from the pole to the equator. It is assumed that meridional flow 
is zero at the pole and equator. The grid is divided so that there are 37 gridpoints 
in the radial direction and 22 gridpoints in the meridional direction, with the radial 
grid size slowiy increasing with radius. The meridional grid is divided so that 
poiints lie equi-distant on either side of 9 = 0 and 9 = 90°, at 9 = -2.25°, 


2.25°, 6.75°, . . ., 87.75°, 92.25°. The algorithm adopted here is the Full-Implicit 
Continuous Eulerian (FICE) scheme described b v Hu a nd_.Wu-Cl.984-). Forrime 
stepping a second-order accurate toward"~difFerencing scheme is used, with the 
step size being of the same order as' given by the Courant condition because the 
magnetic field is calculated explicitly. Smoothing is used when gradients become 
too large, i.e.. at shocks (which do not occur here). At the inner boundary, the 
flow is subsonic and sub-Alfvenic so that two of the six independent variables are 
calculated using compatibility relations (Hu and Wu, 1984). A brief summary of 
the compatibility conditions for the present model is given in the Appendix! along 
with details on how the boundary values and conditions are applied. .We choose 
to specify the radial and meridional magnetic fields, temperature, and density. The 
radial and meridional flow speeds are computed from compatibility relations (i.e-.. 
Equations (A.l) and (A.2)). At the outer boundary, the flow is restricted to being - 
both supersonic and super-Alfv^nic. In this case, all variables at that boundary can . 
be calculated by simple linear extrapolation from the first (or first two) grid points 
inside the boundary. In this study, we did not perform the comparison between 
the present boundary conditions and conventional boundary conditions. However, 
in a recent study by Sun (1991), it was shown that the statement of the boundary 
conditions tn the Appendix eliminates the spurious waves generated by boundary 
disturbances and which can cause numerical instability. _ .... — />. • -■ ' 7 

We start with an essentially initial state and allow the flow to relax in time while 
holding the boundary values'constant. In the present case the initial flow field is 
a polytropic, hydrodynamic solution to the steady-state radial flow equation of 
motion (e.g., Parker, 1963) superimposed on a potential magnetic field. That this 
is neither a self-consistent nor stable solution to the steady-state MHD equations 
is irrelevant since the flow is allowed to evolve in time under the control of the 
equations of motion. TheTpam The main concerns are that the numerical solution 
be stable and of sufficient accuracy to define the physically interesting aspects of 
the solution, and that the relaxation proceed long enough that an acceptably close 
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Fig. 1. Density, temperature, and velocity profiles in radius that were used for the initial (t = 0) 
state in the relaxation. On the left are the proriles for the 3 O^casss and on the right are the 
profiles tor the J = 0.2 case. Note that, except for the velocity scutes, the scales differ between the 
two panels. Because the polytropic index is near unity, temperature changes siowlv with radius and 
the irregularities in tne temperature profiles should be interpreted as noise. 


approximation to tne steady state has been reached. We address these issues briefly 
tn Section 4. 


3. Detailed Results from Four Specific Models 

We report here on four specific models. The results are grouped first according 
to the way in which the physical variables are plotted (i.e., either versus radius 
or versus polar angie) and second according to which of the four examples the 
plot is for. In these four examples, three magnetic field geometries are used: a 
dipole, a quadrupole, and a hexapole; the scalar potentials are therefore proportional 
to Pi (cos 0). P} (cos 0), and P 4 (cos (9), respectively, where P n (cos 9) is the 
Legendre polynomial of degree n. There are two dimensionless free parameters: 
the polytropic index, 7, and the plasma, 3. We use 7 = 1.05 in all cases, 3 = 0.5 
for ail three field geometries, and, in addition, do a dipole calculation for j = 0.2. 
In these case, 3 is evaluated at 1.0 R,t, at the equator, where the field strength is 
1.67 G both for 3 = 0.5 and 3 = 0.2. For the high 3 cases, the base temperature 
and density are 1.8 x 10 6 K and 2.25 x 10 8 cm~^. For the low 3 case, they are 
1.44 x 10 6 K and 5.61 x 10 7 cm" 3 . The three magnetic field geometries naturally 
lead to a single equatorial streamer, a mid-latitude streamer, and both an equatorial 
and a mid-latitude streamer for the dipole quadrupole and hexapole, respectively. 

Results from the four examples will be referred to as follows: 

(a) Dipole, 3 = 0.5: 

(b) Quadrupole, 3 = 0.5« 


OJHG1MAL PAGE \i 
OF POOR QUALITY 


•\ TWO-DIMENSIONAL MHD GLOBAL CORONAL MODEL: STEADY- STATE STREAMERS 5 

(c) Hexapoie, J =:0.5.. 

(d) Dipoie, j = 0.2. 

The initial state temperature, density, and velocity profiles are shown in Figure 1 . 
The temperature curves appear irregular due to the small change in temperature 
over the relatively large radial range - a consequence of the polytropic index being 
near unity. Only three significant figures were retained after the calculation so what 
is seen here is roundoff error in the plotted results rather than in the computed 
results. 

The final, steady-state magnetic field geometries for the four cases are shown in 
Figure 2. Here is seen the weil-known property that the flow is nearly radial beyond 
3-4 R^. The flow is field-aligned everywhere and field lines which cross the outer 
boundary reach to co. The streamers are those volumes which are magnetically 
closed (i.e., the fieid lines return to the surface of the Sun) and it is evident that 
relatively small volumes in the streamers remain magnetically closed in comparison 
to the initial state where all field lines were closed. These closed volumes are 
surrounded by a low density shell but. as will be shown below, the densities in the 
large coronal hole-like open regions are otherwise only slightly lower than in the 
streamers. In each panel of Figure 2. four dashed lines are shown and labelled A, 
B. C. or D. These lines indicate the radial directions used below to plot variables 
versus radius. 

The physical times allowed for the relaxation in these four examples were: 
(a) 22.22 hours for the J = (k5idipoie; (b) 1 6.67 hours for the J = -0 : 5 / 'quadrupole: 
(c) 18.06 hours for the J = 0.5',hexapoie; I'd) 19.44 hours for the J = 0.2 dipole. 
These times are determined by how long it takes for any fluctuation to be advected 
out through the outer boundary of the solution domain. This in turn depends on how 
large the flow speed is and whether the fluctuations represent inward propagating 
waves. In general, the times listed above are the minimum required for a stationary 
fluctuation (i.e., non-propagating in the solar wind frame) to be advected from 
l Rq to 1 5 Rq at a typical flow speed in the open regions. This sometimes leads to 
small residuals in the relaxation near the outer boundary at 15 Rq, but the solutions 
inside 7 R ^ that are shown here are quite steady. This is another point that will be 
reviewed in Section 4. 

Figures 3 and 4 are plots ot density and radial velocity versus radius. The plots 
are made in the directions indicated in Figure 2 so that, for example, in each panel 
of Figure 3 the density is plotted in the four directions A, B, C. and D indicated 
in the corresponding panel of Figure 2. In both of Figures 3 and 4, the four panels 
corresponding to the four panels in Figure 2 are clearly labeled. The densitv profiles 
have been divided by their corresponding initial state (t = 0) profiles from Figure l 
because the density changes by several orders of magnitude between the Sun and 

R/T). The plots here extend only to 7 Rq because there is no new information 
contained outside this radius - the flow is already supersonic and essentially radial. 

Turning briefly to each figure individually, we begin by noting that a density 
enhancement is indicated by values graterthrThnity, and vice versa. The density 
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0.01 f 

1 2 3 4 5 6 7 3 

Solar Radii 



1 2 3 4 5 6 7 8 

Solar Radii 


Fig. 3. Density as a function oradius. Each panel is for the corresponding case in Figure 2cas * 
labelled. The-c^rves are plotted along the directions shown in Figure 2. For example, the four curves 
tor the d = O^exapote labelled (A, B, C. D), are along the four directions shown rn the third panel v. 
of Figure 2 and labelled in the same manner. Each curve has been divided by the initial profile (see 
Figure I). A density enhancement is indicated by values greater than unity, and vice versa. The 
density concentrations in the streamers are clearly visible, generally being on the order of 25% to 
50% above the initial state. 




concentrations in the streamers here are clearly visible, generally beins on the order 
of 25% to 50% above the initial state. The base density for the J — 'O^ases is close 
to that reported by Allen (1955) for the base of the quiet corona and the density 
profile shown here has generally the right behavior for streamers - as shown bv 
curves CJhr cases fa), (c), and fd), and curve D in case (b). Curve D for case (a)A 
the d =v^ 1 5^dipoie. is an example of the density deficit on the flank of a streamer 
that is typical of the results for all the examples. In contrast, the density in the - 
centers of the open regions (curve B in all cases, curve C in case (b), and curve D 
in case (c)) is little different from the initial state, being only slightly smaller. This 
is only surprising when comparison is made to coronal hole observations (Munro 
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Dipole, 3eta=0. 5 . HexaDoie.Beta = 0.5 



12345073 12345673 



12345673 12345673 

Solar Radii Solar Radii 


Fig. 4. Radial velocity as a Junction of radius. Each panel is for the corresponding case in Fieure 2. 
The curves are plotted along the directions shown in Figure 2. as in Figure 3. The velocity inside the 
streamers is seen to be essentially zero. 


and Jackson. 1977) wherein the density was reported to be more than an order of 
magnitude less than in streamers. This difference is a natural consequence of the 
properties of a polytropic model and the choice we have made for the boundary 
conditions on temperature and density — that they be independent of polar ancle. 
The choice leads to both the high density shown here and the low flow speeds 
shown below on open field lines, irrespective of the open streamline geometry. 
To model true coronal hole flow'wiht'a polytropic gas would require at least an 
elevated temperature in the open regions and probably also a lower density at the 
base (Suess et a/., 1977; Suess, 1979). 

The radial velocity is shown in Figure 4, at the positions indicated in Figure 2. 
As described above, and as is generally the case in polytropic models, the flow 
speed in the open regions is similar to the undisturbed initial flow speed shown 
in Figure 1. In the streamer, the flow speed is essentially zero and it is reduced 
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on the adjacent open field lines cue. apparently, to the strongly inclined flow 
direction more than to geometry. The nonzero outflow above streamers te.s., at 
R = 7.14 R.T , at the equator of the 3 = ojjdipoie) refers to the open rield region 
above the streamer’s cusp. The 3 = O.Zaipole is the most extreme example of 
this - and the flow speed is nearly identical to the initial speed everywhere except 
on closed field lines, directly above the center of the streamer, and on the hishiv 
inclined field lines immediately adjacent to the streamer - where the difference is 
still rather small. 

We do not plot the temperature since, due to the polytropic index being 1 .05. 
it varies by only a few percent throughout the computation domain. However, 
this is an effective temperature’ because a polytropic energy equation with a 
polytropic index of 1.05 is equivalent to a large amount of energy being added to 
the flow. Nowhere is the form of this energy specified, nor what the conversion and 
dissipation mechanisms are. However, it has been shown that a polytropic index 
on the order 1 .05 is required to reproduce observations of coronal densities (Suess 
et al., 1977). 

Finally, the magneticaily open regions, although euaivalent to coronal hole 
flows, do not simulate coronai holes because the flow speeds are far too small. 
To obtain reasonable flow speeds in this model it would be necessarv to have 


the temperature vary across the base of the open region - which is well within 
the capability or the model. Such a variation has been shown to reproduce ail 
the known properties of coronal hole flow and lead to accurate simulations of the 
geometry, with the effective temperature being larger in the center of the hole than 
at the edge (.Suess et al., 1977). In contrast to the open regions, the densities in the 
closed regions are similar to observed streamer densities and we feel this model is 
therefore a good approximation to streamer geometry. The temperature must still 
be qualified as an effective temperature, but can be used for diagnostic purpose! in ~ 
combination with planned observations on SoHO/UVCS. 

Some of the results ,cna:be'better viewed and more easily understood when 
plotted versus polar angtelit different heliocentric distances, than versus r a diu s at 
constant polar angles. Such plots are shown for the dnesityT radial velocity, and 
total field strength in Figures 5, 6, and 7, respectivelvT ~~~ 

Figure 5 shows the density drop adjacent 'the streamer. In the panel for the 
3 ^Oojdipole. this drop is quite large, well resolved, and leads into the density 
enhancement inside the equatonal streamer. The only place this does not occur 
is at the base - where the density is held constant. The width of the density 

enhancement in the streamer decreases with height, just as the width of the streamer 

itself decreases with height (e.g.. Figure 2). Essentially the same thin is seen for the 
3 = 0.2 dipole wiht the following quantitative differences: (i) The streamer is much 
higher and wider, (ii) The density depletion on the flanks has a smaller amplitude. 
These differences arejlte primary reason we conclude that solar streamers are better 
described by a 3 =^oJ^>lasma than by a 3 — 0.2 plasma. Qualitatively, a similar 
result is found for the quadrupoie and hexapole. However, it is obvious that the 
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Fig. 5. Density versus polar angle, between the pole (0°) and the equator <90°). Each of the curves ~ 'TTf** c 

is labelled according to the heliocentric distance it refers to/TIuSbthe curves labelled 1.70 R indicate. — — . — 

the density at 1.70 Rq heliocentric radius. The density at thVbase is constant and so the curves there - 

are flat. Above the base, there is a small density enhancement in the streamer (ca. 5% to 50%) and 

a trough in density at the edge of the streamer. In the middle of the open region, the density is very 

close to what it was in the initial state (see also Figure 3). The reason it is not small is that we have 

used constant temperature and density at the base. To produce a true coronal hoie-like profile would 

have required at least an increase in ihe temperature at the base of the open recion (Suess et ai 

1977). 


hexapole is only marginally resolved with the present grid density - there is really 
only one meridional grid point inside the mid-latitude streamer at any given height.* 
The radial velocity in Figure 6 drops precipitously from the magnetically open 
region to the inside of the streamer. Th at th e velocity is not identically zero inside 
the streamer is a result of numerica^iffusirrand is a measure of this numerical 
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artifact in the FICH scheme. For example, at 1.70 R. z . in the J =(05)dipole, the y' 
velocity drops from ca. 60 km s _ v'rs hardly above the noise level in the plots and the 
associated kinetic energy is too smail to affect the dynamics of the solution. Such 
‘slippage’ will, nevertheless, occur in all numerical solutions. At larger heights 
(e.g., 4.90 and 7.14 Rr,) there is small, but finite flow near and in the neutral sheet 
dividing regions of opposite magnetic polarity. This is qualitatively like what is 
observed in the solar wind in the interplanetary medium. The 0 = 0.2 dipole again 
exhibits properties unlike the Sun in the sense that the very low flow speeds inside 
the streamer seem to still exist even at 7.14 Rq - far outside the observed extent 
of closed streamers. 

Figure 7 shows the variation of the total magnetic field strength, {B 2 -f B])'/ 2 , 
across the streamers. The most interesting thing to note in these plots is the enhance- 
ment in total field strength on the flanks of the streamers. This is what ‘confines’ 
the streamers. The field strength for the 3 = 0.2 dipole is seen to vary smoothly, 
with little distinct evidence of the streamer. This is just another indication that the 
presence of the plasma has had little effect on the field geometry in this low-J case. 

4. Accuracy and Stability of Calculations 

This numerical model has been found to be weakly subject to the Courant condition 
on size of time step. Therefore, the size of the time step decreases as the largest 
values of the temperature and magnetic field increase - along with the maximum 
sound and Alfvdn speeds anywhere in the grid. Counteracting this, the higher 
characteristic speeds lead to a somewhat faster relaxation time. However, generally 
shorter time steps are required for smaller /3 calculations. The flow speed also 
plays an important role in determining the relaxation time to a steady state - the 
initial state is a disequilibrium configuration. This imbalance must have time to be 
advected from the base through the outer boundary. The physical time this takes can • 
be estimated by taking a typical (but small) value for the flow speed and caculating 
how long it would take the plasma to flow at this speed from the base to the outer 
boundary. For example, at 150 km s" 1 , to 15 Rq, this takes 18 hours (relaxation 
times we have used here are given in Figure 1). 

A second consideration is gridpoint resolution. The grid used in these examples 
is 4.5 in latitude and about 0.24 R ^ in radius near the base — increasing slowly 
with radius. This is sufficient to adequately resolve the geometry and flow on the 
scale shown in Figure 2. However, if finer scale information is required in, for 
example, the core of the streamers, a denser grid would be required. 

Always a serious consideration in these time-dependent, non-Cartesian MHD 
calculations is the conservation of magnetic flux - that V B = 0 is maintained 
at all times. The condition is maintained here through accurate differencing rather 
than a self-correcting scheme. No anomalous acceleration due to errors in flux 
conservation is apparent in the results. The numerical scheme is pressure-based so 
it is limited by stability to large and moderate j values (e.g., j3 > 0.1) - which 
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Polar Angle (degrees) Polar Angle (degrees) 

Fig. 6. Radial velocity versus polar angle, between the pole and the equator. Each curve, plotted 
for different heliocentric distance, is labelled in the same manner as in Figure 5. The velocity in'... 
the magnetically closed regions is essentially zero. The reason it is not identically zero -is that there - 
is a small amount of numerical diffusion - quite small as indicated by the velocity being less than 
10 km s' inside the 3 ^(^S^iipole streamer at 2.30 


turns out to be the same restriction for maintaining T" B = 0 to the required 
degree. 

Finally, the energy equation 

reduces to v • V ( p / p 1 ) = 0 when a steady state is reached, which means that ( p j p y ) 
is then a streamline constant. This becomes an analytic test of the achievement of a 
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steady-state solution in our case. The boundary values of p and p are the same at all 
latitudes. Therefore, (p/ p 7 ) = 0 has the same value everywhere in the computation 
regime as it has on the boundary if a steady state has been reached. We have checked 
this for the cases shown in Figure 2 and find that for the dipole and quadrupole it is 
constant to within a maximum of 1% and for the hexapole it is constant to within a 
maximum of 4% (average values over the whole grid are less than 1 % in ail cases). 

5. Discussion 

The new feature of this model, with respect to analogous simulations, is the exten- 
sion of the outer boundary to 15 FL?,. This is not a conceptual advance, but this and 
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the stability and ruggedness of the code make it very useful for simulating realistic 
coronal conditions. We present new resuits for quadrupole and hexpoie fields, with 
their accompanying mid-latitude streamers and open magnetic field regions. The 
Alfvdn speed ranged between 800 km s _1 and a few tens of km s -1 . This is lower 
than is believed appropriate for the corona (Suess, 1988), but we expect our model 
will now enable simulations with higher Alfv^n speeds. *»' " 

When comparing our results to those of Steinolfs on, Su ess; and Wu (1982; 
henceforth referred to as SSW), an interesting/ an dimportant d ifference becomes 
apparent. In the present calculation, we have held the density and temperature 
constant at the base, allowing the velocity (and, hence, the mass flux) to ‘float’ 
with time in accordance with the compatbility relations determining the velocity 
from the solution inside the computational domain. In contrast, SSW hold the 
temperature and velocity constant at the base and allow the density to change 
according to the compatibility relationships. SSW determine the location of the 
streamer by locating close field lines and allowing the velocity to decrease to zero 
at the feet of these field lines. A consequence is that inside the streamer, the final 
density is considerably higher than the initial density and this is the primary reason 
for the quantitative differences between their results and ours. 

There is an important consequence of this difference in boundary conditions 
between SSW and the present calculation: the plasma 3 is computed using the 
temperature, density, and magnetic field at the equator and at 1 FL$. This is invariant 
in the present calculation, but in SSW this number is different in the final, steady 
state than at the beginning: there 3 was computed using the initial values. Therefore, 
in SSW in the steady-state solution is actually larger than stated for each example 
they did. Thus, our calculation for a dipole with 0 = 0.2 (case (d)) corresponds to 
cases for 3 < O.LinJSWW We feel that the way we have done the analysis more 
closely corresponds to what occurs and what is physically known for the Sun and •• 
therefore leads to a more precise definition of the problem. So, we conclude that the 
present study has demonstrated a preferable treatment of the boundary conditions 

in comparison to earlier calculations. f I 'C* ^3 

A consequence of the precise examples we have one in cases (a) through (d), ; 
with constant temperature and density, is the flow spee d an TIhigh . "density in the. ' 
magnetically open regions - in comparison to what is believed to be the case 
in solar coronal holes. This is a natural consequence of using a polytropic gas 
in which the flow speed is strongly dependent on base temperature. It also does 
not reflect suggestions from analysis of Skylab data that densities at the base of 
coronal holes may be a factor of two smaller than at the base of streamers (G. Noci, 
private communication). In a continuation of this study, we will produce models 
with varying temperature and density at the base. The variation in temperature 
will, because it is an ‘effective temperature’, reflect a difference in energy balance 
and distribution between the base of coronal holes and streamers instead of a true 
temperature difference. 
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Appendix. The Inner Boundary Conditions According to the Projected 
Normal Characteristic Method; a 2D Case . / > 

The inner boundary conditions are obtained according to the - method of projected 
characteristics (Nakagawa, Hu, and Wu, 1987) with.thpFICE algorithm (Hu and 
Wu, 1984). For the two-dimensional case, the Alfv6nic mode does not exist, there- 
fore, there are six eigenvalues. These six eigenvalues lead to six projected normal 
characteristics and to six compatibiliity equations (see Wang, 1992; for deriva- 
tion). At the inner boundary, since tv > 0 and v T < V„ V f , the characteristics 
dr/df = v r - V s and dr/dt = v r - Vf are towards the lower boundary from interior 
(i.e., outgoing) and need to be considered. There are four incoming characteristics 
(v r , v r + V s , v r + Vf, and one that is degenerate because of the model symmetries), 
so four variables can be specified at the boundary. Two other variables need to be 
calculated from related compatibility equations. We choose the values of B r , Be, 
p, and T to be specified, leaving two quantities (i.e., tv and ve) to be computed 
according to following compatibility equations: 

* 

dvr _ V,B _ + VfC- . : 

dt ~ pV,V f (V} - Vf) ’ . - (A:1) 

dn = vjvA - VJ)B- - Vf{yj - vJ)C - - - . 

dt V s V f {Vj - V})B r Be ’ . ( ) 


with the corresponding variables simplified in two dimensions as follows: 

2 _ .2 _ 

V A = K = — , 

P 

as = ~(RT , 

b 2 = (Bl + Bl) 

P 

Vj = la 2 + b 1 + [(a 2 + 6 2 ) 2 - 4a 2 6 2 ]’^ 2 
l/ J 2 = Ia 2 +6 2 + [(a 2 + 6 2 ) 2 -4a 2 6 2 ] 1 / 2 , 

8 - = P{Vj - Vl)Vf{v r - V f - B r BeV f (vr - v f)^~ 


(A.3 >-; 

(A.4) 

(A.5) 

(A.6) 

(A.7) 
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-< V / - l D - V t)% - BtVfiVr - Vf)2§*~r 

+fj{Vj - Vl)V,V , - B,v,V f aB » _ 
f A ' 1 t ao t de 

-{pa 2 {Vj - V' A 2 ) + B~) 4- BrB 9 V f v e \~- 
- MV? - V?) + ^Ml) ig» _ 


r v ; A/ ' p \rde r dd 
~~~{-P(V? - v?) - B}\ - + B $ ctg$)~ 

- l ^L [p{ y? - K 2 ) - B j ] ctgd - 

T t 

PVfVf , 2 n -2 

— ( v f - V \) + — + P9(Vf - V A )V> . 

C_ = pty- - V'r)V,iv, - V,)2jjz - B r B»V,{v r - V,)^~ 

-<vi - V’) (v, - V',)| + B,vfr, - V,)2£, 

+ri.vl - v/)v,v.i^ - 

r dp r at? 

-(pa J (V A 2 - V-) -B}) + ~- 

- ».(V A J - V?) + M*Z>1 I|£ + Mf»L + 

L p J r<?0 r <90 

+2^W t 2 - Vl) + B|| - + B. ctg«)+ . 

+^W V, 2 - Vl) - B 2 ] ctg * + 


(A.9) 


Since the ideal MHD equations have been used, flow is parallel to the magnetic 
field lines. Thus we determine Be from the relation B r vg = v r Bg\ -- • 
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